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We study the dynamics of the contact-process, one of the simplest nonequihbrium stochastic 
processes, taking place on a scale-free network. We consider the network topology as annealed, i.e. 
all links are rewired at each microscopic time step, so that no dynamical correlation can build up. 
This is a practical implementation of the absence of correlations assumed by mean-field approaches. 
We present a detailed analysis of the contact process in terms of a Langevin equation, including 
explicitly the effects of stochastic fluctuations in the number of particles in finite networks. This 
allows us to determine analytically the survival time for spreading experiments and the density 
of active sites in surviving runs. The fiuctuations in the topological structure induce anomalous 
scaling effects with respect to the system size when the degree distribution has an "hard" upper 
bound. When the upper bound is soft, the presence of outliers with huge connectivity perturbs the 
picture even more, inducing an apparent shift of the critical point. In light of these findings, recent 
theoretical and numerical results in the literature are critically reviewed. 

PACS numbers: 89.75.Hc, 05.70.Jk, 05.10.Gg, 64.60.an 



I. INTRODUCTION 

The study of the effects of an heterogeneous topol- 
ogy on equilibrium and nonequilibrium dynamical pro- 
cesses has lately experienced an active interest from the 
statistical physics community [l|. Indeed, it has been 
observed in recent years that many natural and man- 
made systems are well characterized in terms of complex 
networks or graphs 0, 0|i which vertices represent 
elementary units in the system, while edges stand for 
pairwise interactions between elements. Most real net- 
worked systems can be characterized by a heterogenous 
complex topology, showing remarkable universal features, 
such as the small world property Q and a scale-free con- 
nectivity pattern 0. The small- world property refers to 
the fact that the average distance {£) between any two 
vertices — defined as the smallest number of edges on a 
path between one and the other — is very small, scaling 
logarithmically or even more slowly with the network size 
N This is to be compared to the power-law scaling 
(£) ~ TV^/'' in a d-dimensional lattice. Since the loga- 
rithm grows slower than any power-law function, even if 
d is very large, small-world networks can be thought of as 
highly compact objects of infinite dimensionality. On the 
other hand, scale-free (SF) networks are typically char- 
acterized by a degree distribution P{k), defined as the 
probability that a randomly selected vertex has degree 
k — is connected to k other vertices — that decreases as a 
power-law, 



P{k) 



(1) 



where 7 is a characteristic degree exponent, usually in 
the range 2 < 7 < 3 [3,[1]. 

Dynamical processes taking place on top of complex 
networks arise in a wide variety of scientific and tech- 



nological contexts. For example, we can mention the 
transmission of information packets on the Internet 0], 
the spreading of biological diseases on social networks 
or computer viruses in computer infrastructures [1, 
etc. The interest in the study of these dynamics was 
triggered by the observation that the heterogeneous con- 
nectivity pattern observed in SF networks with diverging 
degree fluctuations can lead to very surprising outcomes, 
such as an extreme weakness in the face of targeted at- 
tacks aimed at destroying the most connected vertices 
ig. ITli . or the ease of propagation of infective agents 
m These properties are due to the critical inter- 

play between topology and dynamics in heterogeneous 
networks and are absent in their homogeneous counter- 
parts. After those initial discoveries, a real avalanche 
of new results have been put forward, including classi- 
cal equilibrium systems [isl . [3, and non-equilibrium 
processes such as epidemic spreading [l^ , reaction- 
diffusion processe s Il6l . [itI [18| and dynamics with ab- 
sorbing states [19I . l2Cll|. For an extensive review of recent 
results we refer the reader to Ref. [ij. 



The analytical approach to the study of dynamical pro- 
cesses on complex networks is dominated by the applica- 
tion of heterogeneous mean-field theory [l|. Heteroge- 
neous mean-field theory (HMF) is based in two basic as- 
sumptions: (i) the homogeneous mixing hypothesis, stat- 
ing that all vertices with the same degree (within the 
same degree class) share the same dynamical properties; 
and (ii) the assumption that fluctuations are not rele- 
vant, and therefore analytical studies can be conducted 
within a deterministic approach. This last fact is in some 
sense natural, since the small- world property implies that 
dynamical fluctuations in a network are so close together 
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that they can be washed away in very few time steps"'^. 
HMF has proved to be extremely useful in providing a 
very accurate description of the behavior of most dy- 
namical processes on complex networks [ij . On the other 
hand, in other instances, such as in the nonequilibrium 
contact process (CP) I24I a debate has arisen about the 
comparison between numerical simulations on SF net- 
works and HMF predictions [H, [H, . Underlying this 
controversy is the fact that while the HMF approach con- 
siders all relevant quantities as deterministic, and hence 
assumes an infinite system size, numerical simulations are 
performed on finite systems and thus are necessarily in- 
fluenced by stochastic fluctuations due to the finite num- 
ber of particles, inparticular close to an absorbing state 
phase transition |22l |. Finite size effects are very strong 
in networks and an appropriate theoretical framework for 
them is necessary to compare simulations with HMF re- 
sults. For this reason, in Ref. [2^ the CP was considered 
on the simplest network substrate, which is an annealed 
network, in which the quenched disorder imposed by the 
actual connections in the network is not considered. In 
this scenario, it was possible to deduce, by means of qual- 
itative arguments, the correct size scaling of the CP in 
this kind of networks, in very good agreement with nu- 
merical simulations. 

In this paper we present a more detailed analysis of the 
CP in SF networks, deriving the corresponding Langevin 
equation describing its dynamics for the case of annealed 
networks. The analysis of this equation allows us to un- 
cover the correct finite size scaling behavior of the CP 
in heterogeneous networks, providing the exact value of 
the critical exponents describing this system. Surpris- 
ingly, the critical behavior of the CP turns out to be 
extremely sensitive to the particular degree cutoff chosen 
for the construction of the network, in agreement with 
previous results obtained from a more phenomenological 
approach [25| . In particular, critical exponents depend 
explicitly on the way the degree cutoff diverges with the 
system size, if it scales sufficiently slowly. If the scaling 
is instead fast, additional complications arise and fluc- 
tuations of the degree distribution strongly perturb the 
picture. 

We have organized our paper as follows. In Sec. [Til we 
describe the main properties of annealed networks, which 
represent the simplest network substrate for a dynami- 
cal process, in which mean-field theory is supposed to 
be exact. We focus, in particular, on the effects of the 
maximum degree allowed on the network cutoff and on 
its fluctuations. Sec. IIIII defines the CP on complex net- 
works, whose mean-field analysis is reviewed in Sec. IIVI 
In Sec. |V] we comment on the different approaches fol- 
lowed in the past to deal with finite size effects on the 



CP in SF networks. The general Langevin theory for the 
CP in networks is presented in Sec. IVI[ while Sec. IVIII 
focuses on the analysis of annealed networks. Sec. IVIIII 
discusses the meaning of finite size effects and finite size 
scaling in heterogeneous networks. In Sec. lIXI wc present 
a digression to the case of annealed networks with out- 
liers, that is, vertices with a degree much larger than the 
average maximum degree expected in the network. Fi- 
nally, we draw our conclusions in Sec. |Xl Some technical 
questions are developed in several Appendices. 



II. ANNEALED SCALE-FREE NETWORKS 

The topological properties of any complex network are 
fully encoded in its adjacency matrix , taking the value 
aij = 1 if there is an edge connecting vertices i and j, 
and zero otherwise. In the so-called quenched networks, 
the values of the adjacency matrix are fixed in time. For 
large quenched networks, a statistical characterization in 
terms of the degree distribution P{k) and the degree cor- 
relations P(fc'|fc), defined as the conditional probability 
that a vertex of degree k is connected to a vertex of de- 
gree k' [2^, [13] , is useful as a compact way to express the 
essential features of the adjacency matrix^. Quenched 
networks are the typical output of most network models, 
such as the configuration model (CM) [2^, [2^, [s^, HH , the 
uncorrelated configuration model (s^ , the class of models 
with hidden variables [ssj , linear preferential attachment 
models 0, [s^l, etc. In this case, each network must be 
considered as a representative of a statistical ensemble 
of random networks, which is characterized by the P{k) 
and P{k'\k) probability distributions. When a dynam- 
ical process takes place on top of such a network, one 
is considering the network as frozen, with respect to the 
characteristic time scale td of the dynamics. In this case, 
in a numerical analysis of a dynamical process, one must 
consider the dynamics over many different quenched net- 
works, all belonging to the network ensemble with the 
same statistically equivalent topological properties, and 
perform an ensemble average to compute the average dy- 
namical quantities. 

In other instances, on the other hand, the very net- 
work is a dynamical object, changing in time over a cer- 
tain time scale r^v. In this case, the correct topologi- 
cal characterization is strictly statistical, given in terms 
of the degree distribution P{k) and the degree correla- 
tions P{k'\k). In the limit tm <C to, that is, when the 
network connections are completely reshuffled between 
any two microscopic steps of the dynamics, while keep- 
ing fixed P(k) and P{k'\k), the resulting networks are 
called annealed (ssl . [Sy, [131 . Apart from the cases where 
they describe the actual evolution of real systems, an- 



At variance with what happens in regular lattices below the crit- 
ical dimension, where in particular, close to a critical point, dy- 
namics is governed by fluctuations [2]| . 



^ A more detailed characterization can be made using higher order 
degree correlations, see Ref. Wh . 
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nealed networks are extremely important from a theo- 
retical point of view, because mean-field predictions for 
dynamical processes on networks are usually obtained in 
this limit, via the so-called annealed network approxima- 
tion In practice one replaces the adjacency matrix 
aij by its ensemble average a(fci, kj), defining the proba- 
bility that two vertices of degree ki and kj are connected. 
This average is given by 



a{k,k') 



NP{k) NP{k') 



k'P{k\k') 

NP{k) 



(2) 



where notation i d k means summation for all vertices 
of degree k. Taking the case of uncorrelated networks, 
withP(fc|A:') = kP{k)/{k) [H, the simple form a(A:,fc') = 
kk'/N{k) results. 

From a numerical point of view, the simulation of dy- 
namics on annealed networks implies the re-generation of 
the whole network every time a microscopic dynamic step 
is performed [s^l • For uncorrelated networks this can be 
efficiently implemented in CP-like dynamics. In this case, 
an annealed network of size N is completely defined by 
its degree sequence {ki, . . . jk^}-, where the degrees ki 
are integer random numbers, extracted according to the 
degree distribution P{k), and restricted between a lower 
bound m and an upper bound M < N. Degree corre- 
lations are given by P{k'\k) = k' P{k') / (k) . Thus, every 
time we need to find a nearest neighbor of a vertex, it is 
selected at random with probability k' P{k')/ (k) among 
the TV vertices present in the network. 

Finite SF networks are additionally characterized by 
another parameter, the degree cutoff kc{N) [111, that is 
the average value of the actual maximum degree fcmax 
in a single realization of the degree sequence: kc{N) = 
(fcmax)- In general, kc is a non decreasing function of the 
network size and, as we shall see below, the CP dynamics 
is very sensitive to its actual size dependence. 

Notice that the value of the cutoff plays a relevant 
role in the determination of degree correlations in finite 
quenched networks [39| . It is known that for the network 
to be closed without degree-degree correlations and no 
multiple edges or self-loops one must impose that degrees 
are smaller than the structural cutoff ~ N^^^ [s^l- In 
uncorrelated annealed networks, however, since they are 
by construction uncorrelated, such a restriction docs not 
apply, and any cutoff is in principle possible. 

Simple considerations based on extreme value theory 
[39| give the probability distribution -Pmax(fcmax) of ob- 
serving a maximum degree fcmax among N degrees inde- 
pendently sampled from a distribution P{k) ^ k^^ and 
bounded by the constraint m < k < M . In the continu- 
ous degree approximation, the distribution of maximum 
degrees takes the form 

P.nax(fc.nax ^ k) ^ N [j l)^-^-^-—-^k~^ (3) 

Using this expression one can compute explicitly the 
value of kc{N), obtaining two different behaviors, de- 



pending on whether M/m is larger or smaller than 
7V1/(t-i), namely 



kc{N) 




M. ^ /vi/(7-i) 

m 



(4) 

where r(z) is the Gamma function [40|. 

For the network to be SF the upper bound of the 
degree distribution must diverge with the system size: 
M ~ N^^^. The parameter a; > is in principle arbi- 
trary, but its value strongly affects the nature of the ac- 
tual maximum of the degree sequence. If M diverges not 
faster than Ari/(7-i) (i. e. w > 7 - 1), then k^ = M is a 
hard cutoff, with no degree larger than fcj,- For w < 7 — 1 
instead, kc ^ 7Vi/(7-i) jg j^g^ ^n average cutoff, but 
(^max^) grows with the network size as Afi+(3-7)/'^ ^ /j2^ 
indicating that fluctuations diverge. 

As a consequence, the maximum degree present in the 
degree sequence has wide fluctuations and outliers, i.e. 
nodes with a degree much larger than kc, may be present 
in the network. It is important to stress that taking 
= 7 — I is very different from setting M ~ 00 from the 
beginning (w = 0) or M = (w = 1), as it is usually 
done in the quenched conflguration model [3l|. In both 
cases the average cutoff kc scales as iVi/(T~i), but for 
a; = 7 — I this is a hard cutoff and fluctuations of the 
value of fcmax around kc are bounded. 

The presence of outliers and large fluctuations in the 
maximum degree has a strong effect on the dynamics 
on annealed networks. In particular, as we will see, the 
relevant quantity characterizing the size effects on the 
dynamics is the second moment of the degree distribution 
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(5) 



The fluctuating nature of this quantity from sample to 
sample can be assessed by looking at its standard devia- 
tion <7g, that can be easily computed, given the uncorre- 
lated nature of the degrees in annealed networks. Thus, 
we have the relative fluctuations 



I 

N 



{k± 

(fc2)2 



- I 



(6) 



Assuming that (fc") ~ (fc„ 



)i+l-7\ 



we have 



^ (fc.nax'"'^) i N'^'^-'^ii-T^), ior COKJ-I 

9" ^ Nkk^J-'^Y ^ \ iv(^-i)/"-\ for > 7 - 1 

(7) 

Thus, fluctuations vanish in the large size limit for to > 
7 — I, while for lo < 7 — f , the fluctuations of g diverge 
as a power law with the network size A^. 

In the rest of the paper, we will mainly discuss the 
simplest case uj > 7 — I, considering often the cases uj = 2 
and w — 7 — 1. The more delicate issue of the effect of 
outliers on the behavior of CP on SF networks will be 
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touched only in Sec. IIXI Notice that the value w = 2 has 
no special meaning here and it is just an example of what 
occurs for > 7 — 1. At odds with the case of quenched 
networks, the structural cutoff kc = N^^'^ does not play 
any role in annealed networks. 



III. THE CONTACT PROCESS ON COMPLEX 
NETWORKS 

We consider the contact process (CP) [12 on hetero- 
geneous networks, which is defined as follows [l^. An 
initial fraction po of vertices is randomly chosen and oc- 
cupied by a particle. Dynamics evolves in continuous 
time by the following stochastic processes: Particles in 
vertices of degree k create offsprings into their nearest 
neighbors at rate A/fc, independently of the degree k' of 
the nearest neighbors. At the same time, particles dis- 
appear at rate p, that, without loss of generality, is set 
to /i = 1. From a computational point of view, the CP 
can be efficiently implemented by means of a sequential 
updating algorithm [H, [13]: At each time step t, a par- 
ticle in a vertex i is chosen at random. With probability 
p = 1/(A -I- 1) the particle disappears. On the other 
hand, with probability 1 — p = A/(A -I- 1), the particle 
may generate an offspring. In this case, a vertex j, near- 
est neighbor of i, is selected at random. If j is empty, a 
new particle is created of it; otherwise, nothing happens. 
In any case, time is updated as < ^ t + [(1 + X)n{t)]~^, 
where n(t) is the number of particles at the beginning of 
the time step. Notice that the factor (1 + A) in the time 
update is due to the fact that each infected particle can 
perform two independent actions, either infect a neigh- 
bor (at rate A) or become healthy again (at rate fi = 1). 
This factor was neglected in previous imp lementations of 
the CP in complex networks [l^, [13, l25l |. The results of 
these works remain, however, unaltered, since the factor 
is irrelevant for steady state properties and amounts only 
to a rescaling for time dependent properties. 

In Euclidean d-dimensional lattices, the CP undergoes 
a noncquilibrium phase transition [22l | between an ab- 
sorbing state, with zero particle density, and an active 
phase, with average constant density of particles, which 
takes place at a critical point Ac. This phase transition 
is characterized in terms of the order parameter p, de- 
fined as the average density of particles in the steady 
state. Defining A = A — Ac, we observe for A < 0, and 
in infinite lattices, an absorbing phase with p = 0. For 
A > 0, on the other hand, the system sets in an active 
phase with a nonzero order parameter, obeying p ~ A^. 
Close to the critical point, the system is also charac- 
terized by diverging correlation length and time scales, 
namely ^ ^ jAj""^ and r ^ |A|~''ii. The critical expo- 
nents /3, i'± and characterize the steady state proper- 
ties of the transition. It is also possible to look at the 
time dependent behavior at the critical point. Thus, for 
example, the particle density is observed to decay in time 
as p{t) ~ . Different quantities can also be defined to 



evaluate the time properties of spreading experiments, in 
which the dynamics evolves starting from a single parti- 
cle. In this case we can define the survival probability, 
S{t), as the probability that the activity lasts longer that 
t, finding at the critical point S{t) ~ . These and 
other critical exponents are not independent, but are re- 
lated by a set of scaling and hyperscaling relations [2^ . 
Thus it is possible to give a full characterization of the 
phase transition of CP in Euclidean lattices using only 
three exponents, that we can take to be (without lack of 
generality) (3, h'± and . Below the critical dimension 
dc = 4, the exponents are nontrivial, and depend explic- 
itly on d. For d > d^ the exponents take the classical 
MF values /? = i^y = 1, i^j. = 1/2. 



IV. HETEROGENEOUS MEAN-FIELD THEORY 
FOR THE CP 

Heterogeneous mean-field theory (HMF) is the basic 
starting point to obtain an analytical understanding of 
the behavior of any dynamical process on a complex net- 
work [l|. In order to take into account the possible fluc- 
tuations induced by the network connectivity, the partial 
densities pk (t) of occupied vertices of degree k d, [4lj are 
considered, from which the total density of particles is ob- 
tained as p{t) = y^,f pk(t)P(k). In the spirit of standard 
mean- field theories |42| , the fact that the quantities pk {t) 
are, in finite networks, of stochastic nature, is neglected. 
Instead, deterministic rate equations are considered, tak- 
ing into account the changes in time of the partial den- 
sities, due to the different steps that the evolution of the 
model can take. 

In the case of the CP, the quantities Pk{t), given by 



Pkit) = 



nk{t) 
NP{k)'' 



(8) 



where ^^(t) is the number of particles in vertices of degree 
fc, can be interpreted equivalcntly as the relative densi- 
ties of particles in vertices of degree k, or the probabilities 
that a given vertex of degree k contains a particle. In a 
step of the CP dynamical evolution, the partial density 
Pk{t) can decrease due to the annihilation of a particle 
in a vertex k (with rate 1), or can increase by the gen- 
eration of an offspring in a vertex k' , nearest neighbor of 
k (with rate A/fc'). Therefore, the rate equations for the 
partial densities in a network characterized by a degree 
distribution P{k) and degree correlations given by the 
conditional probability P{k'\k) can be written as 



dpkjt) 
dt 



~Pk 

(t)-f Afc[l-Pfc(i)]E 



Pik'\k)pk>it) 



(9) 



Given Eq. (O , pfc = is always a solution. The conditions 
for the presence of non-zero steady states can be obtained 
by performing a linear stability analysis (43| . Neglecting 
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higher order terms, Eq. Q becomes 



dt 



}k,k' 



(10) 



exponent can be estimated within HMF, by considering 
the time evolution of the total density at A = 1, namely, 



It is easy to sec that the Jacobian matrix L^ki has a 
unique eigenvector Vk — k and a unique eigenvalue A = 
A — 1. Therefore, a nonzero steady state is only possible 
for A > 0, which translates in a critical threshold for the 
absorbing state phase transition 



Ac = 1, 



(11) 



independent of the degree distribution and the correla- 
tion pattern. 

To get more detailed information on the process, and 
in particular on the shape of the order parameter as a 
function of the rate A, we restrict our attention to uncor- 
related networks. In this case, Eq. ^ reads 



dt 



Imposing the steady state condition, dtPk{i) 
the nonzero solutions 



Pk = 



Xkp/{k) 
1 + Xkp/{k) 



(12) 
0, yields 

(13) 



where pk is now independent of time. By combining 
Eq. with the definition of p, one obtains the self- 
consistent equation for the order parameter p, 



Ap 

(fc) 



fcp(fc) 



Xkp/{k)' 



(14) 



that depends on the full degree distribution. 

In the case of SF networks, for which the degree distri- 
bution in the continuous degree approximation is given 
by P{k) = (7— l)m'^~^k~'^ , with m the minimum degree 
in the network, the solution will depend on the degree 
exponent 7. Substituting the summation by an integral 
in Eq. (|14p . we obtain in the infinite network size limit 
(i.e. when the degree belongs to the range [m, 00]) the 
expression 



p^F 



1,7- 1,7,- 



Xpm 



(15) 



where F[a,b,c,z] is the Gauss hypergeometric fimc- 
tion To evaluate the critical behavior for small p, 

we invert this expression using the asymptotic expansion 
of the hypergeometric function for low densities [40| , ob- 
taining the result p{X) ^ {X - 1)^', with P = 1/(7-2) 
for 2 < 7 < 3 and /3 = 1 for 7 > 3, presenting additional 
logarithmic corrections at 7 = 3. 

Right at the critical point, A = 1, the particle density 
is expected to decay as a power law of time, p{t) ^ t~^ 
[2^ . defining a new, temporal, critical exponent. This 



dpjt) 
dt 



E^(^) 



dpkjt) 
dt 



(16) 



To close this equation, we use a quasi-static approxima- 
tion [l7| . which can be justified in terms of an adiabatic 
approximation for the full Langevin theory for the CP 
(see Sec. lVIBl) . In essence, we consider that, even at the 
critical point, where no steady-state is present, the par- 
tial densities relax to a quasi-stationary state, where they 
take the form given by Eq. (fT3|) . In this case, for SF net- 
works in the continuous degree approximation, Eq. ITF 
will read 



dpjt) 
dt 



-p{t)F 



,7- 



,7- 



(fc) 
p{t)m 



(17) 



Using the asymptotic approximation for the hyper- 
geometric function, valid for low density, we obtain a 
decay exponent in infinite networks given by p{t) ~ t~^ , 
with 9 = (3 for all 7 (logarithmic corrections being again 
present at 7 ~ 3). 



V. FINITE-SIZE SCALING FOR THE CP IN 
COMPLEX NETWORKS 

The exponents obtained within HMF theory in the pre- 
vious Section correspond to the thermodynamic limit of 
a SF network of infinite size. Checking their accuracy 
in numerical simulations becomes thus a nontrivial task, 
particularly close to the critical point, due to the effects 
of finite network sizes. Indeed, because of the small-world 
property, the number of neighbors that can be reached 
starting from a certain node grows exponentially or faster 
with the geodesic distance. This implies that, even for 
large networks, just a few steps are sufficient to probe the 
finiteness of the system. Moreover, in SF networks, lo- 
cal topological properties show very strong fluctuations, 
increasing with the size of the network. 

For general critical phenomena, the theory of finite-size 
scaling (ESS) [11] has successfully overcome this problem 
for processes taking place on regular lattices, allowing 
the detection of the signature of continuous phase tran- 
sitions even in very small systems. For absorbing state 
phase transitions, ESS is based on the observation that, 
even below the critical point, the density of active sites 
in surviving runs ps reaches a quasi-steady state whose 
average is a decreasing function of the system size, and 
that can be expressed as a homogeneous scaling function 
of both the system size and the distance to the critical 
point. In the case of networks, system size is replaced 
by the number of vertices, and the surviving density is 
assumed to fulfill the relation [49, WOl 



p,{A,N)=N-P/''f{AN'^n, 



(18) 
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where f{x) is a scahng function that behaves as f{x) 
for X ^ oo and f{x) ^ const for a; — *■ 0. 
Mean-field theory for homogeneous networks predicts 
the exponents (3=1 and D = 2. For the case of SF, 
Ref. [20[ proposed a phenomenological Langevin equation 
for the particle density, taking the form 



dpjt) 
dt 



Ap{t) - bp{tf ~ dp{t) 



7-1 



(19) 



where ■q{t) is an uncorrelated Gaussian noise. Assum- 
ing a scaling form for the surviving density given by 
Eq. (dH]), and by means of a droplet-excitation argu- 
ment, the authors of [131 found that /3 = 1/(7 — 2) and 
D = (7 — 1)7(7 — 2) for 7 < 3, independent of the network 
cutoff, whenever kc{N) > N^/'' ■2$. 

In Ref. [2^, this issue was pursued by focusing on the 
FSS form of survival probability, which at the critical 
point, and in networks of size N, was assumed to be 



Sit,N)^t-'f{t/t,iN)). 



(20) 



The scaling function f{x) is constant for small values of 
the argument and cutoff exponentially for .t 3> 1. tc{N) 
is a characteristic cutoff time that, according to standard 
mean- field FSS theory should scale as t^iN) ~ N^^^ for 
homogeneous networks [12] . By means of a mapping to a 
biased, one-dimensional random walk, the authors of [25j 
found (5 = 1, while the characteristic time showed the 
form, for heterogeneous networks, tc{N) ^ \J N/g, where 
g is defined in Eq. (O, and is thus dependent on the 
degree cutoff. This surprising result, well confirmed by 
numerical simulations [25l| , is in strong disagreement with 
results of Refs. [10, |23|. in which no cutoff dependence 
was claimed. 

In order to fully ascertain the correct FSS behavior of 
CP in SF networks, we go beyond mean-field and phe- 
nomenological theories and tackle the full problem, tak- 
ing into account its implicit stochastic fiuctuations (par- 
ticularly important in the vicinity of a critical point) by 
means of a Langevin approach. This problem is consid- 
ered in the next Section. 



VI. LANGEVIN APPROACH FOR THE CP ON 
NETWORKS 

A. Generic formalism 

To account for the stochastic fiuctuations of the CP 
close to the critical point, we derive here a Langevin 
equation describing the concentration pk{t) or, alterna- 
tively, the number of active sites of degree fc, nk{t). 
Our derivation follows closely the method developed 
in [13]. We start by deriving exact equations for the 
microscopic dynamics (at the vertex level) of the pro- 
cess. Let (Ti{t) be a random binary variable taking value 
ai{t) = 1 if node i is occupied by a particle at time t 
and (Ti{t) = otherwise. Thus, the state of the process 



at time t is completely determined by the state vector 
Tl{t) = {cri{t), (J2{t), ■ ■ ■ , <^N{t)}- Variables ai{t) can un- 
dergo only two types of transition events: 

1. ai{t) = I ^ ai{t + dt) = 0: Vertex i was occupied 
by a particle at time t, and the particle annihilated 
during the time interval [t, t + dt]. 

2. (Ji{t) = ^ <Ji{t + dt) = 1: Vertex i was empty at 
time t and it received an offspring from an occupied 
nearest neighbor during the time interval [t,t + dt]. 

Assuming that the temporal occurrence of these events 
follows Poisson processes, the previous two events can be 
encoded into a single dynamical equation that describes 
the evolution of ai{t) after an increment of time dt as 

a,{t + dt) = a.,{t)Q{dt) + [1 - a,{t)]T^,{dt), (21) 

where C,i{dt) and rii(dt) are dichotomous random vari- 
ables taking values 

with probability dt 
Udt) - { (22) 

1 with probability 1 — dt 



and 



m{dt) 



1 with probability \dt^^aij(jj{t) — 



with probability 1 



\dty^^ dijCTj {t) — 



(23) 

Eqs. (|22[) and pT|) describe the annihilation of particles, 
while Eq. corresponds to the creation from occupied 
nearest neighbors. 

The set of random variables {C,i{dt)-, i = 1, - ■ ■ , N} are 
statistically independent of each other and of the con- 
jugate random variables {r]i{dt);i = 1, - • • ,N}. On the 
other hand, variables {'iii{dt);i = 1, - • • , iV} are not to- 
tally independent since they may involve common events 
inducing correlations among them. For example, imagine 
two empty vertices, A and B, each of degree 1, connected 
to the same occupied vertex C. Because of the CP dynam- 
ics, during a particle reproduction event at vertex C, the 
particle must choose only one of its neighbors to send the 
offspring. Therefore, if vertex A gets the offspring, vertex 
B cannot receive it and vice-versa, inducing thus corre- 
lations between the random variables rjAidt) and rjB{dt). 
However, it is easy to sec that these correlations are of 
order dt^ and can be then safely neglected. In any case, 
this effect only exists in networks with a quenched topol- 
ogy. In contrast, the annealed network topology changes 
faster than the CP dynamics and, therefore, such corre- 
lations are absent. 

Equations ((2T|) , ((22l) , and ((23)) describe the evolution of 
the state of the system at the most detailed possible level 
of description by specifying the precise state of each and 
every one of the vertices of the network. This description. 
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although exact, is not very useful to derive general prop- 
erties of the system, which are better described by coarse- 
grained quantities. In heterogeneous random networks 
with given degree distribution P{k) and degree-degree 
correlations P{k'\k), the degree of vertices k is the most 
appropriate indicator of the different classes of vertices. 
Therefore, we consider all vertices with the same degree 
to be statistically equivalent. Following these ideas, let 
nk{t) be the number of active vertices of degree k at time 
t, that is, 



Uncorrelated random networks 



(24) 



iGfc 



As we can see, n^-lt) is the sum of a large number of ran- 
dom variables that are nearly statistically independent 
in the quenched version of the network and totally inde- 
pendent in its annealed version. Therefore, by invoking 
the central limit theorem, we expect this variable to fol- 
low a Gaussian distribution and, consequently, to follow 
a Langevin dynamics. To derive the specific form of this 
Langevin equation, we need to calculate the infinitesimal 
moments of nk{t), which can be done using Eqs. (j2ip . 
((22)) . and (|23p . Using the results in Appendix [Xl we can 
finally write the corresponding Langevin equation for the 
CP on annealed networks, namely 



dt 



-nk{t) + A [1 - pk{t)] J2 P{k\k')nk'{t) (25) 
6 it) . lukit) + A [1 - pk{t)] J2 P{k\k')nk>{t), 



where pk{t) = nk{t) / P(k)N is the relative density of ac- 
tive vertices of degree k and {Cfe(t), k = 1, • • • , kc} are 
Gaussian white noises (with zero mean and unit variance) 
uncorrelated among them. Eq. (j25p implicitly assumes 
that nk{t) is a continuous variable. This approximation 
is reasonable as long as nk{t) 3> 1, which is usually the 
case in very large systems, when we consider steady-state 
properties. 

Equation (|25p is one of the main results of this paper 
and is also the starting point for our subsequent anal- 
ysis. As one immediately recognizes, the drift term in 
Eq. corresponds to the standard mean-field approx- 
imation derived in Eq. It is easy to see that the 
potential associated with this drift term has a stable min- 
imum whenever A > Ac = 1 which docs not depend on 
the particular correlation pattern given by P(k\k'). The 
position of this minimum corresponds to the steady so- 
lution in the active phase in the thermodynamic limit. 
The diffusion term, on the other hand, points to a pro- 
cess with multiplicative noise that, as we shall see, has 
important implications when the system is close to its 
critical point in finite size systems. 



Finding solutions of Eq. ([25]) for networks with general 
degree-degree correlations is a rather difficult task. In 
this paper, we focus on the simplest (but instructive) case 
of uncorrelated random networks with a given degree dis- 
tribution P{k). For this class of networks, the transition 
probability takes the simple form P{k\k') ~ kP{k)/{k) 
which allows us to write Eq. ([25)1 as 



dpkjt) 
dt 



~Pkit) + \^[l-pkit)]p{t) 



(26) 



where p{t) = Y.kMt)/N = Efe^(fc)Pfc(0 is the global 
concentration of active nodes at time t and we have di- 
vided Eq. (|25p by the number of vertices of degree fc, 
Nk = NP{k). Analogously, we can write a Langevin 
equation for p{t) as 



dpjt) 
dt 



Pit) A 



(27) 




Pk{t)]p{t) ]^k{t) 



where we have defined A = A — 1 so that the critical point 
corresponds to A = 0. 

Eq. (P7)l is not yet a closed equation for p{t) because 
both the drift and diffusion terms involve the partial den- 
sities Pk{t). To close it, we use an adiabatic approxima- 
tion [43| • From Eq. (|77|) we know that close to the critical 
point, A « 0, p{t) is a slowly varying variable. This is 
due to the fact that the first term in the right hand side 
of Eq. ([?7|) is of order higher than p. On the other hand, 
Pk{t) is a variable that relaxes exponentially fast to its 
quasi-equilibrium state since the lowest order in Eq. (|26[) 
is linear in The adiabatic approximation consists 

in neglecting the term dpk (t) /dt in front of pk (t) and as- 
suming that Pk is a stochastic variable that evolves much 
faster than p{t). Thus, setting dpk{t)/dt = in Eq. pH)) 
and solving for pk{t), wc obtain 



Pk{t) 



\kp{t) 



(k) + \kp{t) 



(28) 



^ It is worth mentioning that this separation of time scales be- 
tween the partial quantities p^, and the global one p has also 



been observed in other dynamics like the A + A 
annihilation process [TtI or the voter models |48|1 



diffusion- 
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The noise term in this equation is subdominant due to 
its dependence on the size of the system. Thus, replacing 
the dominant term in the diffusion one, we finally obtain 



Pk{t) 



\kp{t) 




2\kp{t){k)^ 



(fc) + \kp{t) V Nk [{k) + \kp{t)\ 



■Ut). (29) 



In this way, we obtain an expression for the partial den- 
sities pk as a function of k and p{t) only. Replacing this 
expression in Eq. (j27p and keeping only the first order in 



Nf, , we obtain 



dp(t) 
dt 




kP{k) \kp{t) \ 
(fc) (fc) + \kp{t) j 



(30) 



2\kp{t){ky 



Nk [(fc) + Xkp{t)] 



Notice that the sum of statistically independent Gaus- 
sian white noises is another Gaussian white noise whose 
variance is the sum of the individual variances. Thus, the 
diffusion term in the last equation is, indeed, a Gaussian 
white noise. Therefore, we can finally write 



'^P^'^ -pit)iA^XQ[pit)]) + \l^-^A[pit)m, (31) 



dt 
where 



and 



e[p(i)]^E 



fcP(fc) Xkp{t) 
(fc) (fc) + Xkp{t) 



Mpit)] ^ E 



fcP(fc) (fc-)3 

(fc) m + xkp{tw 



(32) 



(33) 



Eq. ([5T|) is now a closed equation for the total density of 
active vertices p which must be solved with an absorbing 
boundary at p = and a reflecting one at p = 1. As 
we can see from Eq. (|3ip . there is an explicit dependence 
on the size of the network N in the diffusion term of 
the Langevin equation. This size dependence, together 
with the specific functional forms of and A[p] will 
determine the finite size behavior of the system near the 
critical point. 



VII. CP IN ANNEALED SCALE-FREE 
NETWORKS 

In this section we focus on heterogeneous networks 
with a power law degree distribution P(fc) ~ fc"''' with 
fc € [to, M], where 2 < 7 < 3, and M = N^/"^ is the 
degree upper cutoff. In particular, we consider the case 
w > 7 — 1, so that the average maximum of the degree 
distribution kc ^ Ad = N^/'^ is a hard cutoff (Sec. HH. 
The case < 7 — 1 will be considered in Sec. IIXI 



Given the form of the degree distribution it is possible 
to evaluate explicitly the functional form of S[p], that 
determines the dynamical properties of the CP. From its 
definition, Eq. (|32p . it is easy to see that in the limit 
of small density 8[p] has two different functional forms 
depending on whether p is larger or smaller than the 
quantity {k)/Xkc, kc being the network cutoff. Thus we 
have (see Appendix |B|): 




^ < ii; region II 



7-2 



<^ p ^1 region I 



where g = (fc^)/(fc)^ and 

C(7) = mT-2r(^ _ l)r(3 - 7). 



(34) 
(35) 



We denote the regime for small p as region II, and the 
regime for larger p as region I. Analogously, we can eval- 
uate the behavior of the function A as 



A[p] = 



1 



l-^^(7)(^ 



region II 



7-2 



(36) 



region I 



with 



C(7) = -m^"^(7 - 2)(7 - l)r(l + 7)r(l - 7). (37) 

As we can see, the correction term in the region I is 
always very small as compared to 1. Therefore, in the 
rest of the paper we consider that A[p] ~ 1. 

In Fig. [T]the behavior of &[p] is evaluated by numeri- 
cally performing the summation in Eq. ([5^ . The linear 
behavior for small densities is very well obeyed. Instead, 
the scaling of the region I is not cleanly observed even 
for the largest network considered (TV = 10®). This is 
due to the fact that region I is surrounded by two slow 
crossovers, one for p w (fc) /kc (where the transition be- 
tween region I and II takes place) and the other for p = 1 
(where Q becomes independent of p). This has the con- 
sequence that some of the theoretical predictions made 
using the simple approximation given in Eq. (j34p are dif- 
ficult to observe except for extremely large system sizes. 

As a consequence of the form of 0[p], the behavior of 
the system at criticality strongly depends on the type of 
experiment performed to probe the absorbing transition, 
see Fig. [2j Indeed, experiments with stochastic trajecto- 
ries exploring the region p ^ {k)/Xkc feel a drift of the 
form 



Xp_ 

(fc) 



7-2' 



type I drift. (38) 



Instead, any experiment such that trajectories mainly 
stay in the region p ^ {k)/Xkc feels a drift term of the 
form 



*[p] = p(A-AV) 



type II drift. 



(39) 



9 




10 10 10 10 10 10 10 10 10 10 
k p/<k> 



FIG. 1: Numerical evaluation of the function 0[p] as a func- 
tion of kcp/{k) for different network sizes. The degree expo- 
nent is 7 = 2.5 and we use to = j ~ 1. In the type II region, it 
is clearly visible a linear behavior, in agreement with Eq. (|34|) . 
In the type I region, convergence towards the theoretical ex- 
pression given by Eq. p4|l is much slower. 




FIG. 2: Example of two different types of experiment to 
study the CP dynamics at criticality performed in an annealed 
network of size A'^ = 10*, 7 = 2.5, and m = 2. The top 
(black) curve is the density decay starting from a fully active 
network. The bottom (blue) curve corresponds to a spreading 
experiment starting from a single active vertex. The left plot 
corresponds to the hard cutolf M = J\[^^^~'~^') and the right 
plot to M = iV^/^ Grey areas depict the values of p in the 
domain p £ [A*'^^, {k)k~^]. In all cases trajectories are for a 
single run in an instance network. 



The quantity g, that will play a fundamental role in 
the rest of the paper, diverges with the cutoff kc as k'^"^ 
for 7 < 3 (and a; > 7 — 1). Notice that for 7 > 3 the 
leading order is linear in both cases. It is also worth 
stressing that if one lets kc diverge, regime II disappears 
and one is left only with regime I, that coincides with 
what is found using HMF on infinite networks (Sec. IIV[) . 
However, in any finite network, when the density gets 
small it is the drift of type II that rules the dynamics. 



Based on the explicit expression of Q [p] we now provide 
a qualitative and quantitative description of the three 
types of experiment that explore the critical properties 
of the CP dynamics: A) density decay at criticality, B) 
spreading experiments and C) surviving runs. 



A. Density decay at criticality 

Starting from a configuration full of active vertices at 
t = 0, the concentration of active vertices is monitored 
as a function of time until the trajectory is trapped at 
the absorbing boundary. Then an average is performed 
over a large number of different runs up to a time such 
that all runs have survived. In this case, after a ini- 
tial timescale t^, the system first feels the type I drift 
and after a crossover time t* , at very low concentra- 
tions, the type II one. Inserting Eq. (|38p into Eq. (|3ip . 
a pure drift of the type I predicts a behavior pi{t) ~ 
with 9 = 1/(7 — 2). Inserting Eq. ([35]) gives instead 
Pii[t) ~ (ff*)""^- The crossover between the two types of 
behavior occurs for a time t* such that pii{t*) ~ {k)/kc, 
i.e. t* ^ kc/{g{k)) ^ {k)k2~^- A third time scale defines 
the survival time of the different runs that, as we will see 
in the next subsection, scales as tc ~ \/ N/g. 

Fig. [2] shows simulation results for this type of exper- 
iment (top curves) in annealed networks with 7 = 2.5, 
m = 2, for w = 2ora; = 7— 1 for a single run start- 
ing from a fully active network. Different colors (grey 
and white) indicate the different regions depending on 
the shape of the drift term. The first thing to notice is 
that, in the case oi lu = 7 — 1, the region that corre- 
sponds to the type I drift is wider as compared to the 
case oj — 2. Nevertheless, even in this optimal case, we 
do not observe cleanly the two different values of the ex- 
ponent 6. For comparison purposes, in Fig.[2]we also plot 
fimctions t^^ and t^^ that would correspond to the pure 
type I and II behaviors for 7 = 2.5. The exponent 9 ap- 
proaches but docs not reach the theoretical value = 2 
even though simulations are performed in networks of 
size N = 10®. The situation in the case of w = 2 is even 
worse because the crossover happens at shorter times and 
the value 9 = 2 is even more difficult to observe. 

The same is observed in Fig. |21 where we show the 
same as Fig. [2] for networks with 7 = 2.5, w = 2 and 
different network sizes but averaging over 100 runs of the 
process over the same instance network. Additional in- 
formation is provided by the local effective exponent of 
the temporal decay as a function of time (Fig. d]). The 
effective exponent decreases initially quite fast, but even 
for N = 10® the crossover to regime II takes place well 
before the asymptotic value = 1/(7 — 2) = 2 is reached. 
Eventually the effective exponent sets to a constant value, 
that is, quite surprisingly, close to 1.2 instead of the ex- 
pected value 9 = 1. 

The reason why we do not see convincing numerical ev- 
idence of any of the two scaling exponents expected from 
the theory is that the separation of time scales between 
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FIG. 3: Evolution of the concentration p{t) starting from a 
fully infected network for different network sizes, 7 = 2.5, 
and Lo — 2. Results are averaged over 100 realizations in 
a single instance network. Long dashed lines correspond to 
the numerical solution of the set of Eqs. (|12p . using as input 
the empirical degree sequence used in the simulations. The 
very nice agreement between both sets of curves justifies our 
approximation. 




t 



FIG. 4; Local effective exponent as a function of time for 
7 = 2.5, = 2 and various values of A'^. The exponent cor- 
responding to the discrete time ti is given by the slope of the 
line joining p{ti-\) and p{ti+i) in log-log scale. 



very slow convergence of 8[p] to its asymptotic shape 
(Fig. [1]). Despite these difRculties, Fig. 3] suggests that 
by increasing the size of the system wc should eventually 
be able to recover the theoretical exponent 9 = 1/(7 — 2) 

Concerning the exponent corresponding to the type II 
drift, = 1, its evaluation from Fig. [3] is more difficult. 
For instance, in the case w = 7 — 1 the scaling of t* and 
tc is the same and, therefore, the exponent 9 = \ can 
barely be observed. In the case = 2, they scale as t* ~ 
iV(''-2)/2 and ~ iV(''-i)/4. Their ratio then goes as 
t*/tc ^ 7V(t-3)/4 ^;^/g fQj. ^ ^ 2.5) which is a very small 
exponent. The direct consequence is that the evaluation 
of the exponent 9 = 1 from these type of experiments 
is too influenced by the effect of the crossover between 
region I and II. 

To clearly sec the predicted behavior corresponding 
to the type II drift, we perform numerical simulations 
with an initial concentration po well below the critical 
level separating regions I and II. In particular we choose 
Po = {k)/2kc. With this initial conditions the dynam- 
ics is ruled by the type II drift from the very beginning, 
leading to the prediction 

Piiit) = (41) 

In Fig. O wc show simulation results for different network 
sizes as compared to the prediction given by Eq. (|4ip . 
The agreement is very good if we consider that the dashed 
lines in Fig. [5] are generated without fitting any parame- 
ter but using the values of g and po used in the simula- 
tions. 

The conclusion is that observing in simulations the ex- 
ponent 9 = 1/(7—2) predicted by the mean-field the- 
ory in the thermodynamic limit is, although in principle 
possible as a pre-asymptotic regime, too difficult from 
a practical point of view, since one should reach net- 
work sizes that are beyond the capabilities of current 
computers. From Fig. [5] one can estimate that in or- 
der to reach an effective exponent close to 2 a network 
larger than N « 10^^ should be considered. On the other 
hand, the behavior predicted by type II drift spans for a 
shorter time as compared to the type I but is, neverthe- 
less, clearly visible, as shown in Fig. [5] 

B. Spreading experiments 



tx, t* , and tc is too weak. To observe cleanly the two 
different regimes the time scales must be well separated: 
tx ^ tc and this requires very large values of N. 

Two additional elements make the observation of the ex- 
pected scaling of type I, 9 = 1/(7—2), even more difficult. 
First, the time for the onset of scaling (see Appendix [C|) 

^ (7-2)7-ir(3-7)r(7-l)' ^ ' 

can be quite large an it diverges for 7 — > 2. Second, the 



Starting from a single randomly chosen active vertex, 
the activity is followed until it decays into the absorbing 
state and the survival time t is recorded. The survival 
probability S{t), defined as the probability that activity 
lasts longer than t, behaves at the critical point as given 
by Eq. (HnH- 

Fig. [5] shows examples of single realizations of this ex- 
periment (bottom blue curves). In this example, though, 
we have selected realizations that survived a time longer 
than 10^, which roughly corresponds to the value of the 
cutoff time tc for this particular 7 and N. In this way 



11 



1 — I 




FIG. 5; Evolution of the concentration p(t) in region II for 
different network sizes, 7 = 2.5, and uj = 2. The initial con- 
centration is po = {k)/2kc and results are averaged over 500 
realizations. Dashed lines correspond to the theoretical pre- 
diction Eq. (|iT|) . 



we can see the domain of p-space that is visited by the 
trajectories of the experiment. Both for = 7 — 1 and 
Lu = 2, trajectories never reach the white area, where 
the type I drift is dominant. They always remain in the 
domain governed by the type II drift. 

The result of Eq. ([20]) can be derived from the Langcvin 
Eq. pip , using standard techniques of stochastic pro- 
cesses theory |47|. In Appendix [Dl we show that, in the 
limit of an infinite network size, we have 

S{t) = lim S{t, N) = l- e~i/* « (42) 

N^QO t 

that is, we recover and exponent S = I for any degree 
exponent 7. 

The value of this exponent implies that the proba- 
bility density function of survival times in infinite sys- 
tems follows a power law of the form ip(t) ~ and, 
therefore, has diverging fluctuations. However, in finite 
size systems, this distribution has a size-dependent cut- 
off time tc{N), and the divergence of the second mo- 
ment of survival times T2 = {t'^) is then cutoff by tc{N): 
T2 = 2 / tS{t,N)dt ~ 2 f-'t^-^ ~ tciN). A calculation 
of this second moment (Appendix [D| leads to the final 
result 




This expression has an explicit dependence on the size 
of the system but also an implicit one through the size 
dependence of the factor g that, as we have shown before, 
can diverge for 7 < 3 with the system size in arbitrary 
ways. This, indeed, results in an infinite number of ways 
to approach the thermodynamic limit [2^. For 7 > 3 



FIG. 6: Scaling of the survival probability S{t) for 7 = 2.5 
and a; = 7 — 1 in a single instance network. 

instead, 17 is a constant and Eq. (|43p reproduces the well- 
known result of homogeneous MF theory [S^ . 

In Reference [25| it was shown that Eq. (PD|) is obeyed 
for Lu ^ 2 with the scaling of tc{N) given by Eq. (|43p . 
while it is not if no bound is imposed on the degree dis- 
tribution. Fig. [6] shows that the scaling (|20|) holds also 
for a; = 7 — 1 with a hard bound. The violation of the 
scaling occurring when w < 7 — 1 has then to do not with 
the average value of the maximum degree kc ~ Afi/(T~i) 
but with the presence of outliers with exceptionally high 
values of k. 



C. Surviving runs 

In this type of experiment, starting from a given initial 
concentration, only those trajectories that have survived 
for a fixed observation time T > tc are kept and used 
to compute an average concentration of active vertices at 
criticality pg. From a numerical point of view, analogous 
information can be obtained by means of a surviving av- 
erage [2^ . made over the surviving representatives of a 
large number of independent runs. 

The motivation for this type of experiment can be 
traced back to the ESS theory. According to this phe- 
nomenological theory, the concentration of active ver- 
tices in surviving runs satisfies the following scaling rela- 
tion 

p,(A, N) = N-^^/^fiAN^/"). (44) 

For SF networks with 2 < 7 < 3, a phenomenological 
approach in Ref. [l^l predicted (3 — 1/(7 — 2) and D = 
(7— 1)7(7— 2) (see Sec. El). The values of these exponents 
are recovered if one considers the type I drift alone. In 
this section, we show that, in fact, FSS only hold for 
heterogeneous networks in the case of w = 7 — 1. Even 
in this case, because of the slow convergence of the type 
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I drift (e.g. Fig. [T]), the FSS theory presented in [20| 
can only be observed for extremely large systems. This 
goes against the original idea of FSS, which is used to 
recover the critical exponents without the need to reach 
very large systems. 

A critical issue in any FSS theory is the computation of 
the exponent (3 /v. According to Eq. (|44|) . at the critical 
point we expect that the concentration of active vertices 
in surviving runs satisfies 



(45) 



In surviving runs experiments, one selects only those tra- 
jectories that have survived for an arbitrary amount of 
time. Therefore, the probability density function that 
there are n active vertices at time t restricted only to 
surviving runs is 



Ps{n,t\no) 



Pin,t\no) 
S{t\no) ' 



(46) 



where p(n,i|no) is the same probability but measured 
for all trajectories, that is, including those that are ab- 
sorbed at the boundary. Notice that with this definition 
J drips{n,t\nQ) = 1. We are interested in the long time 
limit of this p-d.f, Psin) = limt3>i ps(n, i|no). In this 
limit, the concentration of active vertices at criticality 
for surviving runs is just 



PsiO,N) = 



N 



nps{n)dn. 



(47) 



The probability density function p{n,t\nQ) satisfies a 
Fokker-Planck equation, whose solution allows us to com- 
pute the surviving density (see Appendix IE|) 



p,(0,7V)(x 




1 ( chmN 



r 



Ch){k)N 



7-1' (7-l)fc? 



7-1 



(48) 



where erf(2;) and r(a, z) are the error and incomplete 
Gamma functions, respectively; the first (second) term 
in the right hand side come from type II (type I) drift. 
At this point, the result depends on the particular choice 
of Lo. Suppose first that u; > 7 — 1. In this case, both 
the argument of the error function and the one of the 
incomplete Gamma function diverge as — s- 00. As a 
consequence, the contribution of the type I potential is 
exponentially small and only the first integral contributes 
in the thermodynamic limit, yielding the result 



Ps{0,N) cx 



1 



(49) 



In the case of = 7 — 1, the arguments of both the 
error function and the incomplete gamma function are 
constants in the large size limit. In this case the con- 
tribution of both terms is of the same order in 
Ps{0,N) cx Ar-i/(7-i). Nevertheless, since the effective 
potential (f>{n, N) is a monotonously increasing function 
of n, the contribution in Eq. (|48|) of the type I potential is 
always smaller than that of the type II. The physical pic- 
ture is that trajectories stay most of the time in region I 
except for short excursions to region II that give a small 
contribution that, nevertheless, is of the same order in 
TV. Therefore, we can conclude that the behavior given 
by Eq. ([49]) holds in the whole domain u; E [7 — 1, 2]. In 
terms of to, we can finally write that 



9 



3-7 
2uj 



(50) 



This result implies that the conclusions drawn in 
Refs. [2^ [2^ are essentially incorrect, since the scaling 
of ps{0,N) depends explicitly on the degree cutoff. The 
exponent ratio P/P obtained in [lO, H^l is recovered only 
in the particular case w = 7 — 1. 



VIII. THE MEANING OF FINITE-SIZE 
SCALING 

With all these results at hand, we can now discuss 
which is the role, if any, of FSS theory in the context of 
absorbing phase transitions in SF networks. The aim of 
the FSS ansatz is to connect the behavior of the system in 
the active phase — which is independent of the size of the 
system — for A ^ N-^^"^ and the absorbing one — where 
there is an explicit size dependence — for A <^ N^^"^. 
However, the ability to do so relies upon the "natural" as- 
sumption that the laws ruling the system do not change 
when one performs such transition. In the case of the 
CP dynamics in SF networks, we have shown that in the 
absorbing phase the system is mainly ruled by type II 
drift. However, when A is increased, the concentration 
of active vertices also increases and eventually the system 
starts feeling the type I drift. We are then in a situation 
where there is a change of the underlying laws between 
the active and absorbing phases. Consequently, FSS the- 
ory does not work in this case. Nevertheless, there are 
some subtle details depending on the type of cutoff that 
we discuss next. 
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FIG. 7: Scaling of the density for surviving runs for 7 = 2.2 
and tj = 2. Each point is the result 10^ realizations of the 
stochastic process on each of the 10'^ network realizations. 
The dashed line is an interpolation of the data as a guide to 
the eye. 



In the case of t^; > 7 — 1, the Langevin equation de- 
scribing the dynamics for the concentration p is 



dpjt) 
dt 



Pit) (A - X'gpit)) 



2\p{t) 



N 



m (51) 



that holds if A ^ \g{k)/kc- If we perform the change 
of variables Nf,f = N/g and pef = n/Nef, the previous 
equation becomes 

^ = p^.it) (A - AV./ W) + sj^^m- (52) 

Notice that this equation describes the CP dynamics in 
a homogeneous network of effective size N^f- Therefore, 
Pef{^,Ne.f) must satisfy a FFS with exponents (3=1 
and D = 2, that is 



Pe/(A,iVe/) = 




(53) 



Undoing the change of variables we conclude that 
Ps{A,N) satisfy the anomalous FSS 




for A < 



Xg{k) 



(54) 



T his F SS is anomalous in the sense that when A > 
\/g/N then ps(A, A^) ~ A/5 which depends on the size 
of the system through the factor g. Notice also that for 
7 < 3 the factor g/kc can be reasonably large even for 
large system sizes and, consequently, this anomalous scal- 
ing can be observed in a wide range of values of A. Fig. [7] 
confirms the validity of Eq. ([M]) . 



IX. THE EFFECT OF OUTLIERS 

In the previous sections, we have assumed that the 
maximum allowed degree of the network M scales with 
the system size N as N^/^ , with w > 7 — 1, so that 
the average cutoff in the degree distribution kc is propor- 
tional to M and degrees much larger than k^ are simply 
forbidden. 

When M scales faster than iV^/*^''^^-' instead, the cutoff 
degree kc is not a hard but a "soft" statistical value: 
the maximum degree fcmax in a single realization of the 
network is the result of a random process that yields k^ on 
average but has diverging fluctuations: it is still possible 
to find outlier vertices having degrees much larger that 
kc- To investigate which is the role of outliers in the CP 
dynamics, we introduce a minimal toy network model 
with a hard cutoff kc < iV^/^''^^', just as in the previous 
sections, and then we add a single vertex of degree kout = 
aN, with a G [0, 1]. 

From Eq. (pS]) we see that the concentration in the 
outlier vertex is Pfe„„j ~ 1 provided that p ^ (k) / (XaN), 
where the average degree must be computed including 
the contribution of the outlier. This condition is satisfied 
both in density decay and surviving runs experiments. 
In the case of spreading experiments, the condition is 
satisfied only partially since during the beginning of the 
experiment the density is always of the order p ~ N~^. 
However, in the first two types of experiments, the effect 
of the outlier vertex is that the rest of the vertices "see" 
the outlier always active. Since the outlier holds a macro- 
scopic portion of the edges of the system, all attempts to 
make it active occurring along one of its edges are un- 
successful. The net effect is that the system is shifted 
away from its critical point A = 1 and is effectively in a 
sub-critical state. To quantify this effect and to calcu- 
late the position of the new critical point, we separate in 
Ea. ((32l) the outlier's contribution from that of the rest of 
the vertices. This results in 



0[p(O]= E 



kPjk) \kp{t)/{k) 

(fc) l + Xkp{t)/{k) ^ N{k)'' 



(55) 



that is, the outlier has a constant contribution to the drift 
term whereas the contribution of the rest of the vertices 
goes to zero when p approaches zero. Combining this 
result with Eq. (j3ip we obtain the new critical point 



A° = 



1 



1 - kout/N{k) ■ 



(56) 



To check this result, we perform numerical simulations 
starting from a fully active network in three different sce- 
narios (see Fig. [5]). In the first one, we generate a net- 
work with M = N^/^ and set the dynamics to its critical 
value A = 1. Here, A = 1 is the true critical point and, 
as expected, we find a double power law decay towards 
the absorbing state, as explained in section IVIIl In the 
second scenario, we introduce in the previous network a 
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FIG. 8: Density decay starting from a fully active network for 
the model with and without outlier. Network size is = 10* , 
7 = 2.5, and M = N'^^'^ . Results are averaged over 100 
realizations in a single instance network. 



single vertex of degree kout = N and, again, set A = 1. In 
this case, we observe a clear exponential decay, typical of 
a sub-critical regime. Finally, in the third experiment, we 
keep the network with the outlier but we increase the con- 
trol parameter according to Eq. After this correc- 
tion to the critical point, wc observe again a clear double 
power law decay towards the absorbing state, indicating 
that, indeed, Eq. ((56l) predicts the correct critical point 
A2. 

Along the same lines it is possible to understand also 
surviving runs in the same network with an outlier. The 
equation of motion for the density p at the new critical 
point A = AJ! is 

^ = -(A°W(t) + /^eW, (57) 

where 

9'= k^m/{k?- (58) 

After redefining time in Eq. ([57)1 as t' = (A°)^t, we re- 
cover the same type of Langevin equation as for the case 
> 7 — 1, Eq. (|5ip . but with an effective parameter g' 
given by Eq. ([55]) . From here, we readily obtain the FSS 
form for the average density in surviving experiments, 
namely 

where A' = A — A" . In Fig. [S] the validity of the scaling 
form given by Eq. (j59p is demonstrated. 



FIG. 9: Scaling of the density for surviving runs for 7 = 2.5 
and uj = 2, with at outlier of degree kout = A'^. Each point is 
the result 10** realizations of the stochastic process on a single 
network realization. The dashed line is an interpolation of the 
data as a guide to the eye. 



X. CONCLUSIONS 

In this paper, we have presented a detailed analysis of 
the dynamics of the contact process on annealed scale- 
free networks. Using stochastic differential equations for 
this dynamics, we have clarified the behavior of the model 
close to its critical point and, in particular, its finite size 
scaling. Our results indicate that heterogeneous mean- 
field theory — strictly valid for infinite networks — is prac- 
tically unobservablc for the range of sizes that modern 
computers can reach. The dynamics is instead dominated 
by strong finite size effects that give rise to nontrivial 
anomalous effects. Among them, it is worth to notice 
that the scaling of several relevant quantities (like the 
order parameter close to the critical point or the surviv- 
ing times of the dynamics, etc) depends not only on the 
system size iV, as in regular lattices, but also on the up- 
per cutoff M of the scale-free degree distribution which, 
in general, diverges with the system size as M ~ N^^^ . 
The exponent lo is not fixed by the degree distribution 
alone and, in general, can take different values for differ- 
ent network models or even an arbitrary value that we 
can freely choose in annealed networks. This implies that 
the critical exponents of the dynamics are not universal 
but depend on the arbitrary value of oj. 

Our results allow us to understand the origin of the dis- 
crepancy between the phenomenological finite size scal- 
ing theory pro pos ed in Ref. [20| and the numerical re- 
sults found in [24j25j. Indeed, the Langevin equation 
proposed in Ref. [20|, Eq. ([TO)) , misses the crucial point 
that the coefficient h is not a constant but depends on 
the system size, and grows with it. As a consequence, 
when the system is at its critical point and the concen- 
tration becomes small enough, the term hp^ in Eq. p9p 
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becomes more important than the term dp~'~^ , something 
that would not be possible if b was a constant and 7 < 3. 
This change in the dominating term in Eq. (fT9|) at low 
concentrations invalidates thus the FSS proposed in [20| . 

When w < 7 — 1 an additional interesting complication 
arises: the effective average cutoff of the degree distribu- 
tion kc becomes A^i/(t-i), smaller than M, while its fluc- 
tuations diverge as N grows. This means that, depending 
on the specific realization of the degree sequence, some 
outliers (i.e. nodes with a connectivity much larger than 
the effective average cutoff kc) may appear. We have 
shown that a single outlier connected to a macroscopic 
portion of the system has the effect of introducing an 
apparent shift on the critical point position. The investi- 
gation of the role of outliers (and more in general the role 
of diverging fluctuations in the effective upper cutoff) in 
the contact process and other models is a very interesting 
avenue for further investigations. Notice that this case is 
the relevant one for simulations performed without fixing 
an explicit upper cutoff of the degree distribution, a very 
common habit. 

Last but not least, we would like to stress that the the- 
ory and the simulations presented here give a complete 
understanding of the complex behavior of the contact 
process on annealed scale-free networks. Whether or not 
the same picture also holds for quenched topologies re- 
mains an open question calling for further work. 
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of its infinitesimal moments. The first and second (vari- 
ance) such moments inform us about the expected change 
in the process after an increment of time dt and the vari- 
ance of this expected change. More precisely, if we define 
the variable AX(t) = X{t + dt)-X{t) then the first 
infinitesimal moment is defined as 

^ ^ dt^o dt ^ ' 

Analogously, the infinitesimal variance is defined as 

= Hm ^\^^^mm = ^)^ (^2) 

^ ^ dt^o dt ^ ' 

Functions '^\x\ and D\x\ are called the drift and the dif- 
fusion term, respectively. The Langevin stochastic dif- 
ferential equation can then be written as 



= ^{x(i)\ + ^jwmm. (A3) 



where ^(t) is a Gaussian white noise. 

In our case, we are interested in writing a Langevin 
equation for the coarse-grained quantity nk(t) = 
X^iefc '-'^' (^)- This is convenient for two main reasons: 
First, since it is the sum of almost (or totally) indepen- 
dent random variables, we expect the central limit the- 
orem to hold. This guarantees that the corresponding 
noise in the Langevin equation is Gaussian and white. 
Second, in the thermodynamic limit, nk{t)/N can be 
safely assumed to be a continuous variable and, therefore, 
it is justified the use of (stochastic) differential equations. 

To compute the infinitesimal first moment we write 



APPENDIX A: CALCULATION OF THE 
COARSE-GRAINED LANGEVIN EQUATION 

The derivation of the Langevin equation that describes 
a general stochastic process X{t) involves the evaluation 



(nfe(t + dt)\Y.{t)) = ^^it){C^{dt)) + [1 - a,(t)](ry,(di)) 

(A4) 

where we have made use of Eq. ()21|) . Finally, using the 
probability distributions Eqs. ((^ and (^5]). we are lead 
to 



{nk{t + dt)\Y.{t)) ^nk{t)+dt 



k' iEks-JGk' 

Analogously, we can write an expression for the infinitesimal variance as 



(n|(i + dt)\'E{t)) - {nk{t + dt)\'E(t))^ = dt 



nk 



(A5) 



(A6) 
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To derive the previous equation, we have taken into ac- 
count that, since ai{t) are binary variables taking only 
values or 1, (Tf{t) = ai{t) and that (Ji[t)[l ~ (Ti(i)] = 
0;V/:. Terms of order dt'^ have also been neglected. Un- 
der the annealed approximation, we replace in Eqs. (|A5p 
and (|A6p the adjacency matrix by its average value, 
Eq. (12); which allows us to carry out the sums in 
Eqs. (jASp and (jA6[) and, finally, to obtain the Langevin 
equation Eq. (^5)) . 



APPENDIX B: CALCULATION OF Q\p\ IN 
UNCORRELATED SF NETWORKS 



Let us consider the definition of in Eq. ([5^ . 

namely 



fcP(fc) Afcp/(/c) 
(fc) 1 + Afcp/(fc) ■ 



(Bl) 



The evaluation of this quantity in finite networks depends 
on the value of the particle density p. In particular, if 
p ^ {k)/Xkc, where kc is the network cutoff, then the 
denominator in Eq. (jBl[) can be approximated by unity, 
and we have 



kP{k) Xkp _ (fc2) Xp 



(B2) 



where g = (fc^)/(A;)^. On the other hand, outside this 
region we must keep the full denominator in Eq. (jBl|) . To 
estimate 8 in this case, we perform a continuous degree 
approximation, that is. 



e[p] = 

J ni 



kP{k) Xkp/{k) 
(k) 1 + Xkp/{k) 



1,7- 



Xpm 



1,7- 



(fc) 
Apfcc 



7-2 



{B3) 

where F[a, b, c, z] is the Gauss hypcrgcomctric function. 
Using the asymptotic expansions of the hyp ergeometric 
function for small and large arguments [40|, we can es- 
timate the value of Q in the domain {k)/Xkc ^ p ^ 1. 
Within such domain, the second term in Eq. (jB3p be- 
comes an asymptotically small constant as compared to 
the first term, that yields 



e[p] ~ r(7 - i)r(3 - 7) 



Xpm 

TfcT 



7-2 



(B4) 



APPENDIX C: INITIAL TIME SCALE 



To compute the initial time scale ix needed to reach 
region I starting from an arbitrary initial condition po at 
criticality, we consider Eq. (|3ip with the drift term given 
by Eq. ^ and A = 0, namely 



dp{t) 
dt 



C{l) 



P{t) 



7-1 



(CI) 



The solution of this equation is 
p{t) = 



(fc)7-2 



-1/(7-2) 



(C2) 



The asymptotic state p{t) ^ t ^^^'^ independent of 
the initial condition, is reached for times t, such that 



(fc)7-2 

that is, for t > t^, with 



(7-2)C(7) 2-7 



(C3) 



ix = 



(7-2)C(7) (7-2)^-ir(3-7)r(7-l)' 

(C4) 

where we have used the definition of C(7) in Eq. ([35]) 

and (fc) = (7 - 1)to/(7 - 2). 



APPENDIX D: SURVIVAL PROBABILITY 
EQUATION 

Using standard techniques of stochastic processes the- 
ory, we can obtain the partial differential equation sat- 
isfied by the survival probability of the CP dynamics 
at criticality, starting from an initial concentration pg, 
S{t\po), namely [irf 



dS{t\po) 
dt 



= -Po0[po' 



dS{t\po) , pod^S{t\po) 



dpo 



N dpi 



(Dl) 



This equation is the result of integrating the backwards 
Fokkcr-Planck equation in the domain p G [0, 1] and it 
should be solved with the initial condition S{t = 0\pq) = 
1 and boundary conditions 



S{t\po = 0) and 



dS{t\po) 



dpo 







(D2) 



po = l 



that correspond to an absorbing boundary at p = and a 
reflecting one at p = 1 . The survival probability Eq. ([20|l 
can then be evaluated as 



Sit) = S{t\po = 1/N). 



(D3) 



We first start by evaluating the exponent 5. To this 
end, it is only necessary to solve the problem in the ther- 
modynamic limit N 00. However, the above formula- 
tion is not the most appropriate for this purpose, since 
the solution must be evaluated at po = N^^, that is, a 
value that depends on the size of the system. Therefore, 
we perform the change of variables 



no Npo 



(D4) 



where no is the initial number of active vertices, which 
is eventually set to no = 1 and, therefore, is independent 
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of the system size. Using this new variable, Eq. (|D1 
becomes 



dt 



no 
N 



dS{t\n„) , d^S{t\no) 

no j;-^ . (D5) 



dno 



dnl 



Notice that now the limit N —i- oo can be taken in 
Eq. (|D5p . In this limit, the first term in the right hand 
side of Eq. (jPSP vanishes and the process becomes a 
purely diffusive one with multiplicative noise. The so- 
lution is 



5(i|no) = l-e-"«/*«^. 



(D6) 



Setting finally uq ~ I, leads to the exponent 5 = 1 for 
any 7. 

To evaluate the cutoff tc{N), we compute the second 
moment of the survival times, 12(710), starting from hq 
active sites. However, to compute T2 we first need to 
compute the average surviving time ri(no). It is easy 
to see that Ti(no) = /q S{t\no). Using this result 
in Ea. (|D5p . and assuming that trajectories never feel the 
type I drift, yields the following differential equation for 
Ti{no) 



driQ 



g dTi{no) 

— nn 

N dno 



1 



(D7) 



with boundary conditions ri(0) = and T{{N) = 0. The 
solution of this problem is 



Similarly, the differential equation for T2(no) can also be 
obtained from Eq. (|D5[) as 



(Dll) 



73(710) = 2 / tS{t\no)dt. 
Jo 



This results in the following differential equation (involv- 
ing also Ti{no)) 



d^T2{no) g dT2{no) _ 2Ti{no) 



dn^ 



(D12) 



N dno no 
that satisfies the same boundary conditions as Ti(77o 
The solution of this equation is 



27V pvw 

T2(?^o) = / 

g Jo 



due"" / G(Oe~* dt (D13) 



where 



G{t) = 



due"" 



dq 
— f 

q 



(D14) 



In the limit of large N, this expression can be approxi- 
mated as 




^2(710) =770^/^ J e-* G{t)dt, 



(D15) 



proving then Eq. p3l 



Tl(77o) = \l^ 

9 JO 



due" 



t 



(D8) 



When N is very large, the upper limit in the first integral 
becomes very small. Therefore, we take the limit of the 
integrand when u is close to zero, that is. 



r,(77o)c.J_y^ 



which finally leads to 



du[l+u^^ ][- In7t + 7H ] 

(D9) 



Ti{no) ~ -770 In 



Tto 



27V 



APPENDIX E: PROBABILITY DENSITY 
FUNCTION FOR SURVIVING RUNS 

At the critical point, the probability density p{n,t\no) 
of the number of active vertices at time t given that 
the process had 770 active ones at time t = is ruled 
by a Fokker-Planck equation with a drift term ^'(77) = 
— 776[77/iV] and a diffusion coefficient D{n) = 277, that is. 



d_ 

dn 



779 



I X , / ,| M dp{n,t\no) 
p{n,t\no) [77p(77,i|77o)J = — 



(El) 

(DIO) A direct substitution of Eq. ^ into Eq. (|El|) leads to 



d_ 

dn 



nO 



, ,i . _^ 92 dps{n,t\no) , din [5(^|77o)] 
Ps[n,t\no) + [nps[n,t\no)\ = +Ps(n,t\no) — . 



(E2) 



r 



The density ps(77,t|77o) has, by construction, a well- which is independent of the initial condition. By taking 
defined steady state, that we denote by the limit i ^ 1 in Eq. (jE2p . we obtain 



Psin) = limps{n,t\no), 



(E3) 



l977 



776 



Ps{n) + [nps{n)] = KPs{n), (E4) 

077^ 
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where 



of n. Direct integration of Eq. p4|) gives 



lim-ln[^(t|no)] 
t»i at 



(E5) 



Using the result given in Eq. (jD6[) we conchide that k = 
0, meaning that Ps{n) satisfies the potential solution of 
the Fokker-Planck equation [i^l ■ We can then write that 



p,(0,iV)(xl 



N 



-<ji{n, 



with the effective potential 

0(n,iv) = / e 



dn. 



(E6) 



(E7) 



As in the case of the function 8, the potential (f>{n,N) 
takes a different functional form depending on the value 



2N " 



n < 



Xka 



</)(n,iV) = 



C(7) 
7 



(E8) 

At the critical point, A = 1, we can use this result to 
write 



p,(0,iV)(x - 



{k)N 



N 



(E9) 



where subindices I and II refer to which type of poten- 
tial is dominating the integral. In the limit N ^ 1, we 
can evaluate the contribution of each integral, leading to 
Eq. gll). 
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